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Abstract One route to numerically propagating quantum systems is time dependent density 
functional theory (TDDFT). The application of TDDFT to a particular system’s time evolution is 
predicated on U-representability which we have analyzed in a previous publication. In this work, 
we provide new insights concerning lattice U-representability using an newly developed explicit 
solver for the time-dependent Kohn-Sham potential which contrast with implicit solvers studied in 
the past few years. We present and interpret the force-balance equation central to our numerical 
method, describe details of its implementation, and present illustrative numerical results. A new 
characterization of U-representability for one-electron systems is also included. Taken together, 
the results here open the door to deeper theoretical and numerical investigations of the foundations 
of TDDFT. 


Introduction. — Important classes of time-evolution 
algorithms widely employed by chemists and physicists are 
based on reduced descriptions of the wave function. These 
include methods focused on the two-body reduced density 
matrix and the electron density (the diagonal of the one- 
body reduced density matrix). 

Since all interactions of non-relativistic Hamiltonians 
are between at most two electrons, the IV-electron wave 
function, T, contains more information than necessary. 
For this reason, the two-electron reduced density ma¬ 
trix (2RDM) contains enough information to character¬ 
ize properties of non-relativistic quantum systems [1], 
Unfortunately, one must characterize the set of valid 
2RDMs corresponding to a valid IV-electron wave func¬ 
tion; this is known as the N-representability problem [2]. 
The iV-representability problem was proven to be QMA- 
complete [3] highlighting the theoretical difficulty of 
2RDM methods. Nonetheless, there has been successful 
efforts to perform time evolution using 2RDM methods [4]. 

An even more concise description is afforded by the 
ground state one-electron probability density, 

n t = diag(Tr 2 ...Ar|4't)(4'/|) (1) 

which we will simply refer to as the density. The 
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Hohenberg-Kohn theorems dictate that the ground state 
density, n\ 0 , is sufficient to characterize all properties of 
the quantum system [5]. This provides the basis for den¬ 
sity functional theory (DFT). While theoretically com¬ 
pelling, many functionals to efficiently compute properties 
from the density are unknown. Moreover, the universal 
functional necessary for evaluating the energy is unlikely 
to be determined numerically even to only polynomially 
accuracy in the size of the system. Despite the numer¬ 
ous approximations to the universal functional, compu¬ 
tational complexity arguments [6] showed that obtaining 
the numerically exact functional is intractable even with 
quantum computation. 

The corresponding time-dependent result [7] states, for 
sufficiently well-behaved systems, the potential can be 
computed efficiently with access to a quantum computer. 
Unlike the ground state result, the time-dependent com¬ 
plexity analysis relied on the Kohn-Sham (KS) construc¬ 
tion lying at the heart of nearly all practical schemes for 
DFT. In this letter, we return to the analysis begun in 
our previous work [7] using a combination of theory and 
numerics. 

First let us define the general V -representation problem 
associated with the KS system as the task of constructing 
a model system which has the same expectation values on 
selected observables as a target system. U-representability 
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refers to the existence of solutions to this problem when 
different constraints are placed on the model system. The 
time-dependent subset of ^-representation problems con¬ 
siders as input an initial state and the target trajectory 
of selected observables, and the task is to find the cor¬ 
rect time-dependent fields for a specified control Hamilto¬ 
nian. This general framework is not limited to electronic 
systems as illustrated by a study of this problem in the 
context of spin systems [8]. Here, attention will focus on 
the electronic H-representation problem where the tasks 
is to construct a KS system governed by a time-dependent 
potential V(t) such that the KS density matches the den¬ 
sity of a specified interacting many-electron system at all 
times. 

For fermionic simulations, Ref. [9] was first to give a 
constructive solution to time-dependent V -representation 
problem. This was challenged in Ref. [10] where counter¬ 
examples were presented to this construction. These 
counter-examples were largely addressed by Refs. [11,12] 
through a detailed analysis of densities evolving on lat¬ 
tices. In separate work, an implicit solution using a fixed 
point mapping has been formulated directly in the con¬ 
tinuum limit [13-15]. Here, we will present an explicit 
method based on the algorithm analyzed in [7]. 

The paper begins with the force-balance equation, 
then turn towards the implementation details of the 
solver for the K-representation problem. We give some 
numerical examples before discussing single-electron V- 
representability theorems. Finally, an outlook closes the 
letter. 

The force-balance equation. — The non¬ 
interacting V -representation problem requires that 
the fictitious system’s wave function, |<k t ), evolve such 
that its density expectation value, (d> t |h|<f> t ), matches 
a target evolution, n aim {t). The force balance equation 
determines the required instantaneous potential to cor¬ 
rectly construct the KS system. Note that forces enter at 
second order of evolution as anticipated by Newton’s law: 
F = ma. 

The force-balance equation is easily derived from the 
second derivative of the density following the Heisenberg 
equation [7,12]. If we aim for a target evolution, then 
we should have that d^n aim is equal to i($ t \[H, <9th]|$ t ). 
Expanding this commutator, we have two terms, i[T. dt/h] 
and i\y,dtfi\ which we will physically interpret as well as 
given some guidance on numerical implementation. 

We first discuss the acceleration which the forces must 
cause. The free acceleration, also called the momentum- 
stress tensor [9], Q x = i[f,d t n x ] = — ([T, [T, n x ]]) is in¬ 
dependent of the potential operator. To evaluate the 
free acceleration, we use Q = 23?[diag (Tp^T — p^T 2 )] 
with = a}jCLi as the one-body density matrix. Be¬ 
cause we are considering fermions, cqa.J = Sij — ajoq and 
diCLj = —ajCLi. The forced acceleration is given by the 
difference between the free acceleration and the target ac¬ 
celeration S = d^n aim — Q. The forced acceleration then 


determines the forces required from the potential. 

The forces enter through the term i[V. dth] which can be 
recast into two useful forms; one illustrating a connection 
to forces and the other geared towards determining the 
potential. The first form we examine is 

i[V, d t hj ] = ^2(Vj - V k )T kj (a]a k + 4^). (2) 

k 

This form gives a nice analogy to the real space forces 
as F(x) = —W{x). We note that i(<fr|[t7, d t nj]\$) = 
2J2 k (Vj — V k) 9?[Tjy] has the form of a generalized dis¬ 
crete gradient. Here, the real part of Tf- = Ty (\H|1\H) 
includes both the influence of the probability mass at each 
site as well as the underlying spatial metric. 

The second form is more applicable to numerical sim¬ 
ulation: *(<I>|[V',c? t n. 7 -]|<I>) = (*[«r><9 t n]) 1$)^. To 

introduce a simpler expression, it will be advantageous to 
define M.{A)ij = Aij — Sij (J2 k Ajk)- Then, the force- 
balance equation can be expressed as: 

iM^nr]|*> = E (-2X(^[T $ ])) rs I4 (3) 

S 

For consistency with our publication [7], we define K = 
—2M (5ft[T $ ]) as the force-balance operator. For symmet¬ 
ric matrices, A = A T , A 4(A) will have the constant vec¬ 
tor in its null space. The gauge freedom physically stems 
from the irrelevance of the zero of energy. Since we are 
concerned with time-dependent quantum mechanics, the 
constant potential only imprints an unobservable global 
phase on the wave function. 

V - r e present a lion solver. — In a previous publica¬ 
tion [7], we presented and analyzed an explicit solution 
for the time-dependent potential necessary for TDDFT 
provided with the density time-trace of a K-representable 
system. The algorithm was found to scale polynomially 
in all input parameters except for the K-representability 
parameter [7] which diverges when the interacting system 
no longer has a corresponding KS system. 

The algorithm requires, as inputs, the complete time- 
trace of the density and an consistent initial state $ that 
reproduces the initial density and the initial time deriva¬ 
tive of the density. For numerical implementation, the ki¬ 
netic energy in the lattice basis is also needed. A explicit 
solver based on the Runge-Kutta methods [16] updates 
the wave function based on the KS potential computed 
at each timestep. We discuss some novel aspects of the 
implementation next. 

Preparing initial states. Suppose that an orbital 
should have an initial time derivative given by h a,m . Re¬ 
call [7] that for a one-electron wave function, ip, the den¬ 
sity derivative is given by hj = -iJ2 k T ki({a\a k )ip - 
(4 a »)b). Given that we consider a single electron wave 
function, the 1-RDM elements can be defined as (4 a j) = 
l>i) - 
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To assign the phases, roughly speaking, we must solve 
the equation: A <f> = ^ Ah. This is the content of New¬ 
ton’s method. We will describe the modifications needed 
to handle the gauge degree of freedom after deriving the 
Jacobian, = 0 1 . 

Surprisingly, the Jacobian is also given by (minus) the 
force-balance matrix: 


dhj 

9<t>j 


2 ^ ' Tk i 


5 sin {cj) k - fa) 

d<t>j 


— y J- ki\/^k^i COsfbfe 01) [hjk hij ] 

k 

T ij(( a l a j) + <a]a*» 

-$ij Y T ok(( a { a j) + (a]a k )) 

k 

-Ki, 


(4) 

(5) 

( 6 ) 

(7) 

( 8 ) 


Before applying the Newton method, we must account 
for the gauge corresponding to the global phase of the 
wave function. Other manifestations of this gauge degree 
of freedom are 1) the one-dinrensional null space of the 
Jacobian and 2) the constraint that J2 hj = 0. 

We can fix the phase of one of the M wave function 
components in order to fix the gauge. Suppose the fixed 
phase is the last, then we only update <f) g on the M — f 
remaining components. The fixed gauge Jacobian J g is 
equal to the Jacobian on the first M — 1 components. 

Putting it all together, the Newton rule for updating 
the phase vector </> to a new assignment (p is 

0g = $g-Jg 1 (hg-nf m ) ( 9 ) 

Note that dropping the gauge component to get n g and 
h a g m loses no information due to the constraint that 

E nj = °- 

Numerical results with a straightforward implementa¬ 
tion of Newton’s method works quite well provided that 
the wave function’s initial momentum is somewhat close 
to the target momentum. This is consistent the with ex¬ 
pected performance of Newton’s method in other applica¬ 
tion areas. 

As this paper primarily concerns itself with single¬ 
electron one-dimensional test cases, we only briefly dis¬ 
cuss paths towards adapting the previous method to multi¬ 
electron wave functions. Consider a state with occupied 
orbitals for fi = 1,..., N. The total density derivative is 
merely the sum of each orbital i.e. hj = E^ 7 ^' ■ Thus, 
we can apply the Newton method to the first orbital to 
optimize the phase factors associated with ip\. The or¬ 
thogonality constraints for the remaining N — 1 orbitals is 
then enforced by updating the phases of the other orbitals 
appropriately. 


Numerical solutions to the force-balance equation. To 
handle the inversion in spite of the non-zero vector in 
the kernel of the force-balance matrix, we use the trun¬ 
cated singular value decomposition (tSVD). We also tested 
Tikhonov regularization but we found that the tSVD 
works best in the examples tested. The fixed cut-off used 
for the truncation corresponds to the maximum allowed 
C-representability parameter. 

Note that, at least in the case of a spreading wave func¬ 
tion, the additional vectors in the null space are not in¬ 
dicative that the KS system does not exist. Following the 
theorem in Ref. [12], one may think that it would at least 
correspond to the non-uniqueness of the KS. This is triv¬ 
ially true. 

We can understand this non-uniqueness and remedy 
it easily by considering analogies to the reducibility of 
Markov chains. If a Markov chain is reducible, then there 
is a reordering of the sites such that the Markov matrix 
can be written in blocks as 


P = 


Pi 0 
0 P 2 


( 10 ) 


Suppose P\ and P 2 are both irreducible and aperiodic, 
such that each has a unique fixed point labeled 7Ti and 7T2 
respectively. Now, P has several fixed points i.e. (7Ti, 0) t , 
(0,7t 2 ) t , and (7Ti,7r 2 ) T . The interesting state is, of course, 
(7N,7r 2 ) T while the other two are less interesting. 

In the same way, for densities with disjoint support, 
relabelling sites will also give K a block structure. Then 
the non-trivial solution corresponds to the inhomogeneous 
solutions in each disjoint region. So long as each block 
of K is K-representable, the total system remains V- 
representable with a unique non-trivial global solution. 
When there is a small coupling between two nearly dis¬ 
joint spaces, this must be handled with some care as per¬ 
turbation theory is easily applied to eigenvalues but not 
to eigenspaces. 

We used optimization techniques to tackle this problem 
[17]- Consider the problem of solving Ax = b (with A = 
A T € M(7£) and x, b € TV) on subspace defined by P s.t. 
P 2 = P 


r = |AP;r — 6| 2 (11) 



= x T A\x — x T A\bi — bf A\x + b\ + & 2 

+x T A i 2 A 2 \x - x T Ai 2 b 2 - b 2 A 2 ix (13) 

Setting the derivative of r with respect to the vector x to 
zero implies 

x = (Ai + Ai 2 A 2 i) 1 (A 1 bi + Ai 2 b 2 ). (14) 

In our numerical implementation, we found that this 
expression can be improved by noting that A 2 i is a per¬ 
turbation to the matrix A\. We rearranged the expression 
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to achieve better numerical results: 

x = (Ai + A ± x A 12 A 21 ) 1 (b\ + A 1 1 A 12 b 2 ) (15) 

The cost of the numerical inversion depends on the size 
of Ai since two matrix inversions are needed to evaluate 
Eq. (15) whereas the full inversion only requires one ma¬ 
trix inversion. In numerical experiments, we found that 
when dim[Ai] is less than 60% of the full space, Eq. (15) 
is faster. 

Additional improvements were made by employing a 
double truncation technique whereby one domain is de¬ 
fined by region where the density is non-trivial and a sec¬ 
ond domain is defined as the region where n is non-trivial. 
In the examples where the wave function is spreading, the 
region corresponding to h is larger. If we let P' be the 
larger of the two domains and P the smaller, the proce¬ 
dure employed solves for the KS potential within region 
P' and then truncates the potential to region P. 

Another issue that affects numerical performance is the 
choice of the gauge. The two choices we considered were 
1) fixing the gauge such that the mid-point of the domain 
defines zero potential, 2) fixing the gauge as the mean of 
the potential. The first choice is natural in many of the 
examples since the potential is known to be zero at the 
mid-point. This is not generically the case so the second 
choice offers a sensible alternative. By choosing the gauge 
based on the average value of the potential, the norm of the 
Hamiltonian takes its minimal value over all gauge choices. 
This improves numerical stability when the timesteps are 
based on the norm of the Hamiltonian but we also found 
that it can cause erratic jumps in the KS potential due to 
numerical noise at the boundaries. 

Numerical examples. — In this study, we tested our 
solver on some simple numerical examples to validate the 
solver: ballistic spreading, a superposition of particle-in- 
a-box states and the “discovery” of a constant potential. 

The first example, ballistic spreading of the wave func¬ 
tion, served as a simple test of the proposed truncation 
procedure and provided insights into the best methods for 
handling spreading. It was used to compare several differ¬ 
ent truncation procedures and served as validation for the 
procedure presented here. The numerical noise led to a 
potential norm on the order of 10 -8 as the wave function 
spread ballistically. The second example merely confirmed 
the results of [11] and also served as a test case for the pro¬ 
cedure when the density is evolving in a confined space. 

The last test case considered is a density which evolves 
from an initial Gaussian distribution with unit variance 
centered at the origin and zero momentum in a fixed po¬ 
tential defined by 

( -4 \x — 5.5| < 1.5 

V{x) = l -4 |5.5 — x\ < 1.5 (16) 

[ 0 otherwise 

When the propagation starts, the density does not know 
about the potential wells to the right and left, but as it 
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Fig. 1: As the density spreads, it discovers the potential and 
this is reflected in the solutions to the F-representation prob¬ 
lem depicted here. The numerical errors at the boundary of 
the wave packet’s extent are all less than unity and the aver¬ 
age maximum error at the boundary is 0.3195. 

evolves the second derivative and the density itself, learn 
about the potential wells. The numerical procedure re¬ 
mains stable throughout as depicted in Fig. 1. The discov¬ 
ery of the potential is clearly illustrated as the density be¬ 
gins to evolves into the lower potential regions. The worst 
numerical error encountered during the numerical proce¬ 
dure occurs when the density reaches the outer edges of 
the potential wells as shown in Fig. 2. This can be under¬ 
stood as the solver’s inability to decide what the potential 
should be beyond the area that the density has seen. 

Here, it is interesting to note that the spreading of the 
density only occurs at second order. From d/ri = i[H,h] 

d t fij = —i'^ j T kj (a' j a k - a' k a,j) (17) 

k 

Suppose that Idt) = d'lfl) = Y^Gk\K) with K = 

nf=i has no support on site j. Then [ aj , 'F]|fi) = 0 

since 4* has no support on j. Therefore, (ST|'F+a.^.a.ytF|T2) 
is zero as well. Hence, 

(T| d t hj\4/) = -i^2T kj ((0|T t aJa fc ’F|f2) - h.c.) = 0 

k 

(18) 

Thus, density spreads to new regions only at second order 
in time. 

V- r opreselit a bi 1 i ty. — We now make a few com¬ 
ments on the existence of solutions to the F-representation 
problem. A previous theorem [12] showed that ground 
states many-body interacting systems are always V- 
representable in the neighborhood of the initial time. 
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Fig. 2: The difficulty of solving the E-representation problem 
is illustrated by the potential and density at time 1.0671 as 
plotted here. The red dot-dashed line indicates the cutoff used 
for the density and the heavy blue line represents the potential. 

They showed that the matrix K has only one zero eigenval¬ 
ues and is positive definite in the space of inhomogeneous 
potentials. It should be noted that the theorem does not 
characterize the E-representability parameter thus numer¬ 
ical stablility is not ensured. 

Here we use simpler arguments to provide additional 
characterizations of the spectrum of K in the single¬ 
electron case. This coincides with the previous theorem 
for the ground state but generalizes to all eigenstates. 

Theorem: Given non-degenerate ip such that Hip = 
Ekip for H = T + E e M(7£), K(ip) has k — 1 negative 
eigenvalues and M — k + 2 positive eigenvalues. 

Proof: Assuming that ip is an eigenstate of H = T+ V 
with eigenvalue A, then H^ = T — (A1 — E) = T — 
D has ip in its null space. Rearranging, (T + D)ip = 0 
implies TjfcUfc = Djtp-j- Before using the definition of 
the force-balance equation, it is important to note that 
eigenvectors of symmetric matrices are real. Hence, for a 
single-particle in the eigenstate ip: 

Ktj = 2Tijipiipj + 25ijipj (y, Tjktp^j (19) 

= - 2 Tijipiipj - 25ijDjipj ( 20 ) 

= yX S miA)2{Tmn ~ S m nD n )(S nj 1pj) ( 21 ) 

mn 

( 22 ) 

Since ip has full support, S m j = S m jipj is non-singular and 
the number of (+/0/—) eigenvalues are the same for H^ 
and K{ip) by Sylvester’s theorem [18]. □ 

According to numerical tests, if ip does not have full 
support then K will have an additional vector in the null 


space. This is consistent with the theorems from Ref. [12]. 
The interacting extension of the present theorem does not 
seem to hold although we found that the many-body non¬ 
interacting ground state gives rise to I\ with the same 
inertia as consistent with the previous findings [12]. 

Outlook. — Next steps for the long term project be¬ 
gun here are the study of interacting electronic examples, 
comparisons and combinations with the implicit fix-point 
methods [13-15], and designing algorithms for building 
multi-particle initial states. Theoretical questions to be 
tackled include understanding E-representability of open 
system evolutions, continuum limits and intersections with 
quantum computing. Previously, we have shown that the 
E-representation problem can be solved efficiently using a 
quantum computer but it still remains an open question if 
all quantum computations remains efficiently simulatable 
with TDDFT given access to efficient solutions to the V- 
representation problem. 

* * * 
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